## Replication File for Reneging on Alliances: Experimental Evidence Final 
rm(list = ls())
library(rstatix)
library(foreign)
library(tidyverse)
library(ggplot2)
library(MASS)
library(stargazer)
library(gmodels)
library(Rmisc)
library(xtable)
library(cowplot)
library(gridExtra)
library(cobalt)
library(dplyr)

## read the rawdata

rawdata<- read_csv("/Users/victorxu/Desktop/side project/military alliance /research and politics/data replication/Final/raw data copy.csv")

## set the dataset *d1* for analysis in main context (time between 7 minutes to 60 minutes)

d1<- rawdata %>%
  filter(Duration..in.seconds. > 420 & Duration..in.seconds. < 3600) %>%
  glimpse()

## write.csv(d1,"/Users/victorxu/Desktop/side project/military alliance /research and politics/data replication/data for reweighting.csv")

## recode the variable

## recode approval rate as binary variable (our major dependent variable)

d1$biatti<- d1$X7attitude
d1$biatti[d1$X7attitude>=5 & d1$X7attitude <=7]<- 1 ## approve 
d1$biatti[d1$X7attitude>=1 & d1$X7attitude<=4]<- 0 ##disapprove
table(d1$biatti)


## replicate Figure 1: Plot of Treatment Effects (95 percent confidence interval)

summary_table <- d1 %>% group_by(treat) %>% get_summary_stats(biatti, type = "mean_ci") 
newtable<- summary_table [c(6,1,2,5,3,4),]
newtable$name <- c("Control (C)",
                   "Control (C)+ large casualties",
                   "C+ low prospect of victory",
                   "C+ economic justification",
                   "C+ UN mediation",
                   "C+ economic sanctions"
                   
)
newtable$groupno<- c(1,2,3,4,5,6)
newtable
newtable$effect <- ifelse((newtable$name  ==  "Control (C)"), 
                          (newtable$mean - newtable$mean [newtable$name  ==  "Control (C)"]),
                          (newtable$mean - newtable$mean [newtable$name  ==  "Control (C)"]))
newtable



Approve_CI95_3_a <- ggplot(subset(newtable, newtable$name != "Control (C)"), aes(x =name, y = effect*100)) +
  geom_hline(yintercept = 0, linetype = "longdash") +
  geom_errorbar(width = .1, aes(ymin = (effect-ci)*100, ymax = (effect+ci)*100), position = position_dodge(0.25)) +
  geom_point(size = 4, position = position_dodge(0.25)) +
  scale_y_continuous(breaks = seq(-10,20,5), limits = c(-10,20))+ ## change to percentage +
  
  ylab("The effects between control group and treatment groups (% point)\n") + xlab("") + theme_bw() +
  scale_color_grey(start = 0.65, end = 0) +
  scale_shape_manual(values = c(20, 18)) +
  theme(axis.text = element_text(colour = "black", size = 11),
        title = element_text(colour = "black", size = 13, face = "bold"),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        legend.justification = c(1, 1), legend.position = c(1, 1), 
        legend.title = element_blank(), legend.background = element_blank())

Approve_CI95_3_a


Approve_CI95_3_a+scale_x_discrete(name ="", 
                                  limits=c("Control (C)+ large casualties","C+ low prospect of victory","C+ economic justification","C+ economic sanctions","C+ UN mediation"))


ggsave("/Users/victorxu/Desktop/side project/military alliance /research and politics/data replication/figure1.jpeg",
       width = 297, 
       height = 180, 
       units = "mm",
       dpi=300)

### Replicate the simple t.test in the manuscript and regression models in Table A1-A4 in the Appendix 

## recode the treatment

d1$tr1<- NA
d1$tr1[d1$treat==6]<- 0 ## control 
d1$tr1[d1$treat==1]<- 1 ## large casualties
table(d1$tr1)
table(d1$treat)

d1$tr2<-NA
d1$tr2[d1$treat==6]<- 0 ## control 
d1$tr2[d1$treat==2]<- 1 ## low prospect of victory 
table(d1$tr2)


d1$tr3<-NA
d1$tr3[d1$treat==6]<- 0 ## control 
d1$tr3[d1$treat==5]<- 1 ## economic impact
table(d1$tr3)

d1$tr4<-NA
d1$tr4[d1$treat==6]<- 0 ## control 
d1$tr4[d1$treat==3]<- 1 ## UN mediation
table(d1$tr4)


d1$tr5<-NA
d1$tr5[d1$treat==6]<- 0 ## control 
d1$tr5[d1$treat==4]<- 1 ## economic sanctions
table(d1$tr5)

## treatment 1 large casualties and control group 

t.test(biatti~tr1,data=d1) ## p=0.013
t.test(X7attitude ~tr1,data=d1) ## p=0.005
t.test(repu1~tr1,data = d1) ## p=0.44
t.test(repu2~tr1,data = d1) ## p=0.45
## regression models 
m1.1<- glm(biatti~tr1, data=d1, family = "binomial")
summary(m1.1)
m1.2<- glm(biatti~tr1+gender+race+age+inc+educ+pid+employ, data=d1,family = "binomial" )
summary(m1.2)
m1.3<- polr(factor(X7attitude)~tr1, data=d1, Hess = TRUE)
summary(m1.3)
m1.4<- polr(factor(X7attitude)~tr1+gender+race+age+inc+educ+pid+employ, data=d1, Hess = TRUE)
summary(m1.4)
m1.5<- polr(factor(repu1)~tr1, data=d1, Hess = TRUE)
summary(m1.5)
m1.6<- polr(factor(repu1)~tr1+gender+race+age+inc+educ+pid+employ, data=d1, Hess = TRUE)
summary(m1.6)
m1.7<- polr(factor(repu2)~tr1, data=d1, Hess = TRUE)
summary(m1.7)
m1.8<- polr(factor(repu2)~tr1+gender+race+age+inc+educ+pid+employ, data=d1, Hess = TRUE)
summary(m1.8)

## treatment 2 low prospect of victory and control 

t.test(biatti~tr2,data=d1) ## p=0.07
t.test(X7attitude ~tr2,data=d1) ## p=0.015
t.test(repu1~tr2,data = d1) ## p=0.102
t.test(repu2~tr2,data = d1) ## p=0.318

## regression models 

m2.1<- glm(biatti~tr2, data=d1, family = "binomial")
summary(m2.1)
m2.2<- glm(biatti~tr2+gender+race+age+inc+educ+pid+employ, data=d1,family = "binomial" )
summary(m2.2)
m2.3<- polr(factor(X7attitude)~tr2, data=d1, Hess = TRUE)
summary(m2.3)
m2.4<- polr(factor(X7attitude)~tr2+gender+race+age+inc+educ+pid+employ, data=d1, Hess = TRUE)
summary(m2.4)
m2.5<- polr(factor(repu1)~tr2, data=d1, Hess = TRUE)
summary(m2.5)
m2.6<- polr(factor(repu1)~tr2+gender+race+age+inc+educ+pid+employ, data=d1, Hess = TRUE)
summary(m2.6)
m2.7<- polr(factor(repu2)~tr2, data=d1, Hess = TRUE)
summary(m2.7)
m2.8<- polr(factor(repu2)~tr2+gender+race+age+inc+educ+pid+employ, data=d1, Hess = TRUE)
summary(m2.8)

## treatment 3 economic impact and control

t.test(biatti~tr3,data=d1) ## p=0.44
t.test(X7attitude ~tr3,data=d1) ## p=0.96
t.test(repu1~tr3,data = d1) ## p=0.91
t.test(repu2~tr3,data = d1) ## p=0.43

## regression models 
m3.1<- glm(biatti~tr3, data=d1, family = "binomial")
summary(m3.1)
m3.2<- glm(biatti~tr3+gender+race+age+inc+educ+pid+employ, data=d1,family = "binomial" )
summary(m3.2)
m3.3<- polr(factor(X7attitude)~tr3, data=d1, Hess = TRUE)
summary(m3.3)
m3.4<- polr(factor(X7attitude)~tr3+gender+race+age+inc+educ+pid+employ, data=d1, Hess = TRUE)
summary(m3.4)
m3.5<- polr(factor(repu1)~tr3, data=d1, Hess = TRUE)
summary(m3.5)
m3.6<- polr(factor(repu1)~tr3+gender+race+age+inc+educ+pid+employ, data=d1, Hess = TRUE)
summary(m3.6)
m3.7<- polr(factor(repu2)~tr3, data=d1, Hess = TRUE)
summary(m3.7)
m3.8<- polr(factor(repu2)~tr3+gender+race+age+inc+educ+pid+employ, data=d1, Hess = TRUE)
summary(m3.8)


## treatment 4 UN mediation and control 

t.test(biatti~tr4,data=d1) ## p=0.1975
t.test(X7attitude ~tr4,data=d1) ## p=0.1089
t.test(repu1~tr4,data=d1) ## p=0.016
t.test(repu2~tr4,data=d1) ## p=0.1225

## regression 
m4.1<- glm(biatti~tr4, data=d1, family = "binomial")
summary(m4.1)
m4.2<- glm(biatti~tr4+gender+race+age+inc+educ+pid+employ, data=d1,family = "binomial" )
summary(m4.2)
m4.3<- polr(factor(X7attitude)~tr4, data=d1, Hess = TRUE)
summary(m4.3)
m4.4<- polr(factor(X7attitude)~tr4+gender+race+age+inc+educ+pid+employ, data=d1, Hess = TRUE)
summary(m4.4)
m4.5<- polr(factor(repu1)~tr4, data=d1, Hess = TRUE)
summary(m4.5)
m4.6<- polr(factor(repu1)~tr4+gender+race+age+inc+educ+pid+employ, data=d1, Hess = TRUE)
summary(m4.6)
m4.7<- polr(factor(repu2)~tr4, data=d1, Hess = TRUE)
summary(m4.7)
m4.8<- polr(factor(repu2)~tr4+gender+race+age+inc+educ+pid+employ, data=d1, Hess = TRUE)
summary(m4.8)

## treatment 5 economic sanctions and control 

t.test(biatti~tr5,data=d1) ## p=0.1812
t.test(X7attitude ~tr5,data=d1) ## p=0.1899
t.test(repu1~tr5,data = d1) ## p=0.2434
t.test(repu2~tr5,data = d1) ## p=0.4547

## regression 
m5.1<- glm(biatti~tr5, data=d1, family = "binomial")
summary(m5.1)
m5.2<- glm(biatti~tr5+gender+race+age+inc+educ+pid+employ, data=d1,family = "binomial" )
summary(m5.2)
m5.3<- polr(factor(X7attitude)~tr5, data=d1, Hess = TRUE)
summary(m5.3)
m5.4<- polr(factor(X7attitude)~tr5+gender+race+age+inc+educ+pid+employ, data=d1, Hess = TRUE)
summary(m5.4)
m5.5<- polr(factor(repu1)~tr5, data=d1, Hess = TRUE)
summary(m5.5)
m5.6<- polr(factor(repu1)~tr5+gender+race+age+inc+educ+pid+employ, data=d1, Hess = TRUE)
summary(m5.6)
m5.7<- polr(factor(repu2)~tr5, data=d1, Hess = TRUE)
summary(m5.7)
m5.8<- polr(factor(repu2)~tr5+gender+race+age+inc+educ+pid+employ, data=d1, Hess = TRUE)
summary(m5.8)

## Replicate Table A1

stargazer(m1.1, m2.1,m3.1,m4.1, m5.1,m1.2,m2.2,m3.2,m4.2,m5.2, type = "html", 
          title            = "Logit Estimates: Determinants of Public Approval for Non-military Intervention",
          covariate.labels = c("Large Casualties", 
                               "Low Prospect of Victory",
                               "UN Mediation",
                               "Economic Sanctions",
                               "Negative Effect of Economy",
                               "Gender","Race","Age","Income","Education","Party Identity",
                               "Employment"),
          dep.var.caption  = "Approval Rate",
          dep.var.labels   = "Approval Rate",
          column.labels   = c("Bivariate Models",
                              "Full Models"),
          column.separate = c(6, 2),
          omit.stat=c("LL","aic"), no.space=TRUE,
          align = TRUE,
          style = "qje",
          notes = ("Two-tail significance levels; Clustered standard errors"),
          out="/Users/victorxu/Desktop/side project/military alliance /research and politics/data replication/TableA1.html")

## Replicate Table A2

stargazer(m1.3, m2.3,m3.3,m4.3, m5.3,m1.4,m2.4,m3.4,m4.4,m5.4, type = "html", 
          title            = "Ordered-Logit Estimates: Determinants of Public Approval for Non-military Intervention",
          covariate.labels = c("Large Casualties", 
                               "Low Prospect of Victory",
                               "UN Mediation",
                               "Economic Sanctions",
                               "Negative Effect of Economy",
                               "Gender","Race","Age","Income","Education","Party Identity",
                               "Employment"),
          dep.var.caption  = "Approval Rate (7 Point Scale)",
          dep.var.labels   = "Approval Rate (7 Point Scale)",
          column.labels   = c("Bivariate Models",
                              "Full Models"),
          column.separate = c(6, 2),
          omit.stat=c("LL","aic"), no.space=TRUE,
          align = TRUE,
          style = "qje",
          notes = ("Two-tail significance levels; Clustered standard errors"),
          out="/Users/victorxu/Desktop/side project/military alliance /research and politics/data replication/TableA2.html")

## Replicate Table A3

stargazer(m1.5, m2.5,m3.5,m4.5, m5.5,m1.6,m2.6,m3.6,m4.6,m5.6, type = "html", 
          title            = "",
          covariate.labels = c("Large Casualties", 
                               "Low Prospect of Victory",
                               "UN Mediation",
                               "Economic Sanctions",
                               "Negative Effect of Economy",
                               "Gender","Race","Age","Income","Education","Party Identity",
                               "Employment"),
          dep.var.caption  = "Perceived Impact on the US Future Military Alliance Commitments (5 Point Scale)",
          dep.var.labels   = "Perceived Impact on the US Future Military Alliance Commitments (5 Point Scale)",
          column.separate = c(6, 2),
          omit.stat=c("LL","aic"), no.space=TRUE,
          align = TRUE,
          style = "qje",
          notes = ("Two-tail significance levels; Clustered standard errors"),
          out="/Users/victorxu/Desktop/side project/military alliance /research and politics/data replication/TableA3.html")
## Replicate Table A4
stargazer(m1.7, m2.7,m3.7,m4.7, m5.7,m1.8,m2.8,m3.8,m4.8,m5.8, type = "html", 
          title            = "",
          covariate.labels = c("Large Casualties", 
                               "Low Prospect of Victory",
                               "UN Mediation",
                               "Economic Sanctions",
                               "Negative Effect of Economy",
                               "Gender","Race","Age","Income","Education","Party Identity",
                               "Employment"),
          dep.var.caption  = "Perceived Impact on the US Future Non-Military Commitments (5 Point Scale)",
          dep.var.labels   = "Perceived Impact on the US Future Non-Military Commitments (5 Point Scale)",
          column.separate = c(6, 2),
          omit.stat=c("LL","aic"), no.space=TRUE,
          align = TRUE,
          style = "qje",
          notes = ("Two-tail significance levels; Clustered standard errors"),
          out="/Users/victorxu/Desktop/side project/military alliance /research and politics/data replication/TableA4.html")


## Replicate Figure 1 "Univariate Balance on Pre-Treatment Covariates" in the Appendix

library(gplots)    # version 3.1.1
library(extrafont) # version 0.17
## recode
df<- d1 %>% glimpse()
## female 
df$female <- NA
df$female [df$gender==2]<-1 ## female
df$female [df$gender==1]<-0 ## male 
## 
par(mfrow = c(2,3))
par(mar = c(1.8,1.8,1.8,1.8))
## white 
df$white <- NA 
df$white [df$race==1]<-1 ## white
df$white [df$race>=2 & df$race<=7]<-0 ## non-white  
table(df$white)
## college degree 
df$col<- NA
df$col[df$educ<=4]<-0 ## lower than 4-year college
df$col[df$educ>=5]<-1 ## college or above

## employment 
df$emp<- NA
df$emp [df$employ==2]<- 0 ## non-employment
df$emp [df$employ==1]<- 1 ## employment

## republican 
df$republic<- NA
df$republic [df$pid_r==1|df$pid_r==2|df$pid_n==1] <- 1 ## republican 
df$republic [df$pid_d==1|df$pid_d==2|df$pid_n==2|df$pid_n==3]<-0 ## non=republican
## martial 
df$ married<- NA 
df$married [df$marital==1]<-1 ## married 
df$married [df$marital>=3]<-0 ## non-married
table(df$fp1)
# Gender
plot(density(df[df$treat==1,]$female, na.rm=T, bw=0.1), xlab="Gender (0 = Male, 1 = Female)", ylim=c(0,3.5), xlim=c(-0.4,1.4), main="Gender", lty=1, cex=1.2)
par(new=TRUE)
plot(density(df[df$treat==2,]$female, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.5), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=2)
par(new=TRUE)
plot(density(df[df$treat==3,]$female, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.5), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=3)
par(new=TRUE)
plot(density(df[df$treat==4,]$female, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.5), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=4)
par(new=TRUE)
plot(density(df[df$treat==5,]$female, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.5), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=5)
par(new=TRUE)
plot(density(df[df$treat==6,]$female, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.5), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=6)
par(new=TRUE)
par(new=TRUE, family="Times")
legend("topleft", legend=c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5", "Group 6"), lty=c(1,2,3,4,5,6,7,8), cex=1, bty="n")


# Race
plot(density(df[df$treat==1,]$white, na.rm=T, bw=0.1), xlab="White (0 = No, 1 = Yes)", ylim=c(0,3.2), xlim=c(-0.4,1.4), main="Race", lty=1, cex=1.2)
par(new=TRUE)
plot(density(df[df$treat==2,]$white, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.2), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=2)
par(new=TRUE)
plot(density(df[df$treat==3,]$white, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.2), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=3)
par(new=TRUE)
plot(density(df[df$treat==4,]$white, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.2), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=4)
par(new=TRUE)
plot(density(df[df$treat==5,]$white, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.2), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=5)
par(new=TRUE)
plot(density(df[df$treat==6,]$white, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.2), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=6)
par(new=TRUE, family="Times")
legend("topleft", legend=c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5", "Group 6"), lty=c(1,2,3,4,5,6), cex=1, bty="n")

# Education
plot(density(df[df$treat==1,]$col, na.rm=T, bw=0.1), xlab="College degree (0 = No, 1 = Yes)", ylim=c(0,3.7), xlim=c(-0.4,1.4), main="College degree", lty=1, cex=1.2)
par(new=TRUE)
plot(density(df[df$treat==2,]$col, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.7), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=2)
par(new=TRUE)
plot(density(df[df$treat==3,]$col, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.7), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=3)
par(new=TRUE)
plot(density(df[df$treat==4,]$col, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.7), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=4)
par(new=TRUE)
plot(density(df[df$treat==5,]$col, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.7), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=5)
par(new=TRUE)
plot(density(df[df$treat==6,]$col, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.7), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=6)
par(new=TRUE, family="Times")
legend("topleft", legend=c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5", "Group 6"), lty=c(1,2,3,4,5,6), cex=1, bty="n")

# Party ID
plot(density(df[df$treat==1,]$republic, na.rm=T, bw=0.1), xlab="Republican (0 = No, 1 = Yes)", ylim=c(0,3), xlim=c(-0.4,1.4), main="Party affiliation", lty=1, cex=1.2)
par(new=TRUE)
plot(density(df[df$treat==2,]$republic, na.rm=T, bw=0.1), xlab="", ylim=c(0,3), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=2)
par(new=TRUE)
plot(density(df[df$treat==3,]$republic, na.rm=T, bw=0.1), xlab="", ylim=c(0,3), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=3)
par(new=TRUE)
plot(density(df[df$treat==4,]$republic, na.rm=T, bw=0.1), xlab="", ylim=c(0,3), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=4)
par(new=TRUE)
plot(density(df[df$treat==5,]$republic, na.rm=T, bw=0.1), xlab="", ylim=c(0,3), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=5)
par(new=TRUE)
plot(density(df[df$treat==6,]$republic, na.rm=T, bw=0.1), xlab="", ylim=c(0,3), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=6)
par(new=TRUE, family="Times")
legend("topright", legend=c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5", "Group 6"), lty=c(1,2,3,4,5,6), cex=1, bty="n")

# Marital status
plot(density(df[df$treat==1,]$married, na.rm=T, bw=0.1), xlab="Married (0 = No, 1 = Yes)", ylim=c(0,3.3), xlim=c(-0.4,1.4), main="Marital status", lty=1, cex=1.2)
par(new=TRUE)
plot(density(df[df$treat==2,]$married, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.3), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=2)
par(new=TRUE)
plot(density(df[df$treat==3,]$married, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.3), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=3)
par(new=TRUE)
plot(density(df[df$treat==4,]$married, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.3), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=4)
par(new=TRUE)
plot(density(df[df$treat==5,]$married, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.3), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=5)
par(new=TRUE)
plot(density(df[df$treat==6,]$married, na.rm=T, bw=0.1), xlab="", ylim=c(0,3.3), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=6)
par(new=TRUE, family="Times")
legend("topright", legend=c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5", "Group 6"), lty=c(1,2,3,4,5,6), cex=1, bty="n")

# Employment status
plot(density(df[df$treat==1,]$emp, na.rm=T, bw=0.1), xlab="Employed (0 = No, 1 = Yes)", ylim=c(0,4), xlim=c(-0.4,1.4), main="Employment status", lty=1, cex=1.2)
par(new=TRUE)
plot(density(df[df$treat==2,]$emp, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=2)
par(new=TRUE)
plot(density(df[df$treat==3,]$emp, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=3)
par(new=TRUE)
plot(density(df[df$treat==4,]$emp, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=4)
par(new=TRUE)
plot(density(df[df$treat==5,]$emp, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=5)
par(new=TRUE)
plot(density(df[df$treat==6,]$emp, na.rm=T, bw=0.1), xlab="", ylim=c(0,4), xlim=c(-0.4,1.4), main="", axes=F, ylab="", lty=6)
par(new=TRUE, family="Times")
legend("topright", legend=c("Group 1", "Group 2", "Group 3", "Group 4", "Group 5", "Group 6"), lty=c(1,2,3,4,5,6), cex=1, bty="n")

dev.off()